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Abstract 


The  inverse  scattering  problem  for  an  acoustic  medium  is  formulated  by 
using  the  variable  background  Born  approximation.  A  constant  density  acoustic 
medium  is  probed  by  a  wide-band  point  source,  and  the  scattered  field  is 
observed  along  a  curved  receiver  array  located  outside  the  region  where  the 
medium  velocity  is  different  from  the  assumed  background  velocity  function. 

The  solution  that  we  propose  relies  on  the  introduction  of  a  backpropagated 
field.  This  field  is  obtained  by  using  a  finite-difference  scheme  backwards 
in  time  to  backpropagate  into  the  medium  the  scattered  field  observed  along 
the  receiver  array.  The  backpropagated  field  is  imaged  at  the  source  travel 
times,  giving  an  image  of  the  same  type  as  obtained  by  reverse-time 
finite-difference  migration  techniques.  The  greidient  of  this  image  is  then 
taken  along  rays  linking  the  source  to  points  in  the  medium,  and  after 
scaling,  this  gives  the  reconstructed  potential.  To  relate  the  reconstructed 
potential  to  the  true  scattering  potential,  we  use  high  frequency  asymptotics 
and  an  additional  approximation  Introduced  by  Beylkin  [1].  These 
approximations  reduce  the  validity  of  our  reconstruction  procedure  to  the  high 
wavenumber  region.  With  these  approximations,  it  is  shown  that  at  a  given 
point,  the  reconstructed  potential  corresponds  to  the  convolution  of  the  true 
potential  with  a  weighting  function  obtained  by  partially  reconstructing  an 
impulse  from  its  projections  inside  a  cone.  The  angular  range  of  this  cone  is 
totally  determined  by  the  geometry  of  the  receiver  array,  and  by  the  relative 
location  of  the  source  with  respect  to  the  point  that  we  consider.  In  the 
special  case  when  the  receiver  array  surrounds  the  domain  where  the  scattering 
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potential  is  located,  we  find  that  within  the  Born  approximation,  the 
reconstructed  potential  recovers  exactly  the  high  wavenumber  part  of  the 
Fourier  transform  of  the  true  potential.  It  is  expected  that  for  a  wide  class 
of  problems,  the  reconstruction  technique  described  in  this  paper  will  be 
computationally  more  efficient  that  the  generalized  Radon  transform  (GRT) 
inversion  method  proposed  by  Beylkin,  Miller  and  Oristaglio  [1]  -  [3]. 
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1 .  Introduction 

The  multidimensional  inverse  scattering  problem  for  acoustic  media  has  a 
wide  range  of  applications  in  areas  such  as  exploration  geophysics  [4], 
ultrasonic  imaging  [5],  and  nondestructive  testing,  among  others.  In  this 
problem,  the  medium  of  interest  is  probed  by  sources  located  outside  the 
medium  and  the  scattered  field  is  recorded  at  various  locations.  The 
objective  is  to  reconstruct  the  velocity  and  density  of  the  medium  as 
functions  of  position.  In  this  paper,  we  will  assume  for  simplicity  that  the 
density  of  the  medium  is  constant.  Since  the  relation  between  the 
observations  and  the  velocity  function  is  nonlinear,  the  solution  of  the 
inverse  problem  must  also  be  nonlinear.  For  one-dimensional  (1-D)  media,  i.e. 
for  media  whose  velocity  function  varies  only  with  one  space  dimension,  a 
number  of  exact  inverse  scattering  procedures  which  rely  either  on  integral  or 
on  differential  equations  have  been  proposed  over  the  years.  A  discussion  of 
these  methods  can  be  found  in  [6].  However  for  multidimensional  media,  exact 
inverse  scattering  methods  are  still  at  an  early  stage  of  development,  and  the 
methods  which  have  been  developed  until  now  require  experiment  geometries 
which  are  rather  unrealistic  from  a  practical  point  of  view. 

This  has  motivated  researchers  to  develop  approximate  direct  inversion 
methods,  which  solve  exactly  a  linearized  form  of  the  multidimensional  inverse 
scattering  problem.  The  Born  and  Rytov  approximations  [7]  are  the  most 
commonly  employed  of  these  linearization  techniques.  Whether  one  should  use 
one  of  these  approximations  instead  of  the  other  depends  on  the  nature  of  the 
scatterers  and  on  the  experiment  geometry.  In  this  paper,  we  shall  consider 
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the  Born  approximation.  In  this  approximation,  the  object  profile  that  we 
want  to  reconstruct  is  viewed  as  a  small  perturbation  about  an  assumed 
background  velocity  model,  and  the  scattered  field  is  expressed  linearly  in 
terms  of  this  perturbation.  Physically,  the  Born  approximation  takes  into 
account  only  the  singly  scattered  waves;  multiple  scattered  waves  due  to  the 
velocity  perturbations  are  neglected.  Note  however  that  the  multiples  due  to 
the  background  model  are  included  in  the  scattered  field. 

During  the  past  few  years,  a  number  of  solutions  of  the  multidimensional 
Born  inversion  problem  have  been  proposed  for  various  observation  geometries 
and  background  velocity  models.  The  majority  of  these  solutions  assume  a 
homogeneous  background  model.  The  zero-offset  reflection  geometry  consisting 
of  coincident  sources  and  receivers  was  considered  by  Cohen  and  Bleistein  [8] 
for  a  line  aperture  in  two  dimensions,  and  by  Norton  and  Linzer  [5]  for  plane, 
cylindrical  and  spherical  apertures  in  three  dimensions.  More  recently, 
Fawcett  [9]  formulated  the  zero-offset  Born  inversion  problem  as  a  tomographic 
problem,  where  the  objective  is  to  reconstruct  a  function  from  its  projections 
along  circles  or  spheres.  The  inversion  method  proposed  by  Cohen  and 
Bleistein  was  also  extended  to  the  case  of  a  stratified  (1-D)  background  model 
in  [10]  and  [11].  Raz  [12]  and  Clayton  and  Stolt  [13]  considered  the  same 
experiment  geometry  as  Cohen  and  Bleistein,  but  with  unstacked  data.  In  this 
geometry,  some  source  and  receiver  arrays  are  located  on  the  surface  of  the 
earth,  and  for  each  source  the  scattered  field  is  recorded  at  all  receivers 
rather  than  only  at  the  coincident  receiver.  They  showed  that  both  the 
density  and  bulk  modulus  of  the  acoustic  medium  can  be  recovered  with  this 
experiment.  Another  interesting  feature  of  Clayton  and  Stolt’s  paper  is  that 
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it  assumes  a  variable  background  model,  which  is  then  taken  into  accoxint  by 
extrapolating  the  observed  scattered  field  downwards  into  the  medium  of 
interest. 

The  solution  of  the  multidimensional  Born  inversion  problem  which  was 
proposed  by  Esmersoy  and  his  colleagues  [14]-[16]  relies  on  a  similar 
backpropagation  principle.  However,  they  considered  a  different  scattering 
experiment,  where  the  medium  is  probed  by  a  single  wide-band  point  or 
plane-wave  source,  and  where  the  scattered  field  is  measured  along  a  curved 
receiver  array  located  outside  the  region  where  the  medium  velocities  differ 
from  the  background  model.  In  this  approach,  the  scattering  potential  is 
reconstructed  by  propagating  the  observed  scattered  field  backwards  in  time 
with  a  finite-difference  scheme,  imaging  this  field  either  at  zero  time  [15] 
or  at  the  source  travel  times  [14],  [16],  and  then  filtering  the  resulting 
image.  Depending  on  whether  the  receiver  array  surrounds  the  scattering 
object  or  not,  the  reconstructed  potential  is  an  exact  or  approximate  solution 
of  the  linearized  inverse  scattering  problem.  However,  an  important 
limitation  of  the  results  described  in  [14]-[16]  is  that  they  were  restricted 
to  the  homogeneous  background  case. 

The  first  complete  solution  of  the  variable  background  Born  inversion 
problem  was  presented  by  Beylkin,  Miller  and  Oristaglio  [1]  -  [3],  who 
formulated  this  problem  as  a  generalized  Radon  transform  (CRT)  inversion 
problem,  where  the  objective  is  to  reconstruct  a  fxinction  from  its  projections 
along  a  set  of  curves  whose  geometry  depends  on  the  experiment  and  on  the 
background  model.  The  solution  obtained  by  Beylkin  and  his  colleagues  is 
expressed  in  terms  of  a  weighted  backprojection  operator  summing  the 
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contributions  of  the  curves  passing  through  a  given  point.  This  solution  is 
quite  general  since  it  applies  both  to  the  case  when  the  medium  is  probed  by  a 
single  point  source,  and  to  the  zero-offset  £ind  finite-offset  geometries. 

Note  however  that  the  inverse  GRT  relies  on  high  frequency  asymptotics,  and  on 
an  additional  approximation  which  together  have  the  effect  that  the 
reconstructed  potential  recovers  only  the  rapid  space  variations  of  the  true 
potential,  or  equivalently  the  high  wavenumber  part  of  its  Fourier  transform. 
Thus,  the  inverse  GRT  solves  the  variable  background  Born  Inversion  problem 
only  in  a  limited  sense.  Another  limitation  of  the. inverse  GRT  is  that  to 
obtain  the  weights  appearing  in  the  backprojection  operator,  one  must  use  a 
ray  tracing  algorithm  for  every  point  along  the  receiver  array  where 
measurements  are  taken.  The  amount  of  computations  required  by  this  method  is 
therefore  quite  large. 

The  objective  of  this  paper  is  to  extend  the  backpropagated  field 
approach  of  [14]-[16]  to  the  variable  background  case  for  the  experiment  where 
the  medium  is  probed  by  a  single  wide-band  point  source.  The  computational 
procedure  that  we  use  to  obtain  the  reconstructed  potential  is  quite  different 
from  the  inverse  GRT,  but  the  domain  of  validity  of  our  reconstruction  method 
is  the  same  as  that  of  the  inverse  GRT.  In  other  words,  the  reconstructed 
potential  recovers  only  the  rapid  space  variations  of  the  true  potential,  or 
equivalently  the  high  wavenumber  part  of  its  Fourier  transform.  The  first 
step  in  our  approach  is  to  filter  the  observed  time  traces  and  to  use  them  as 
source  wavelets  for  doublet  sources  located  along  the  receiver  array.  For 
this  source  configuration,  a  finite  difference  scheme  is  used  backwards  in 
time  to  compute  the  backpropagated  field  inside  the  medium.  This  field  is 
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then  imaged  at  the  source  travel  times,  and  the  reconstructed  potential  is 
obtained  by  taking  the  gradient  of  this  image  along  rays  linking  the  source  to 
points  in  the  medium,  and  by  scaling  the  resulting  expression.  To  relate  the 
reconstructed  potential  to  the  true  potential,  we  use  high  frequency 
asymptotics  as  well  as  an  additional  approximation  introduced  by  [1].  These 
approximations  reduce  the  domain  of  validity  of  our  inversion  method  to  the 
high  wavenumber  region.  In  this  context,  it  is  shown  that  at  a  fixed  point, 
the  reconstructed  potential  is  the  convolution  of  the  true  potential  with  a 
weighting  function  obtained  by  partially  reconstructing  an  impulse  from  its 
projection  inside  a  cone.  The  angular  range  of  this  cone  is  totally 
determined  by  the  geometry  of  the  receiver  array  and  by  the  location  of  the 
source  with  respect  to  the  point  that  we  consider.  A  consequence  of  this 
representation  is  that  in  the  special  case  when  the  receiver  array  provides  a 
total  coverage  of  the  scattering  region,  and  within  the  domain  of  the  validity 
of  the  Born  approximation,  our  reconstruction  procedure  recovers  exactly  the 
high  wavenumber  part  of  the  Fourier  treinsform  of  the  scattering  potential. 

The  step  of  our  reconstruction  method  which  is  the  most  demanding  from  a 
computational  point  of  view  is  the  computation  of  the  backpropagated  field 
with  a  finite-difference  scheme.  We  expect  that  for  a  significant  number  of 
problems,  this  computation  will  be  easier  to  perform  than  the  ray  tracing 
scheme  which  is  required  for  every  receiver  by  the  inverse  GRT. 

This  paper  is  organized  as  follows.  In  section  2,  the  velocity  inversion 
problem  is  formulated  within  the  Born  approximation.  The  definition  and 
implementation  of  the  backpropagated  field  are  discussed  in  section  3.  Our 
reconstruction  method  is  described  in  section  4,  and  a  representation  theorem 
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is  obtained  for  the  weighting  function  relating  the  reconstructed  potential  to 
the  true  potential.  This  representation  relies  on  approximations  which  reduce 
the  domain  of  validity  of  our  reconstruction  procedure  to  the  high  wavenumber 
region.  The  representation  that  we  obtain  shows  that  at  a  given  point,  the 
weighting  function  is  the  partial  reconstruction  of  an  impulse  from  its 
projections  over  a  cone.  The  construction  of  this  cone  is  examined  in  section 
5,  and  is  shown  to  depend  only  on  the  experiment  geometry.  Section  6  contains 
some  conclusions  and  some  thoughts  about  further  extensions  of  our  results. 
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2.  Problem  Description 

Consider  the  scattering  experiment  described  in  Fig.  1.  A  constant 
density  2-D  acoustic  medium  is  probed  by  an  impulsive  point  source  located  at 
Tj  in  the  x-y  plane.  Since  this  2-D  medium  is  in  fact  a  3-D  medium  whose 
velocity  function  c{x)  does  not  vary  with  the  third  dimension  z,  the  2-D  point 
source  is  in  fact  implemented  by  a  line  source  parallel  to  the  z  axis.  The 
scattered  field  is  observed  along  a  curved  receiver  array  F,  and  in  the 
following  it  is  assumed  that  F  is  a  smooth  curve  parameterized  by  the  arc 
length  s,  sel. 

Note  that  since  we  assume  that  the  2-D  medium  is  probed  by  a  line  source, 
the  experiment  geometry  described  above  is  not  totally  realistic.  In 
practice,  it  is  much  cheaper  to  use  a  point  source  to  probe  the  medium.  In 
this  case  the  problem  becomes  a  2  1/2-D  problem  in  the  sense  that  although  the 
velocity  function  c(x)  varies  only  with  two  space  dimensions,  the  waves 
propagate  in  three  dimensions.  The  geometrical  spreading  associated  to  the 
three-dimensional  propagation  of  the  waves  needs  to  be  taken  explicitly  into 
account  in  the  inversion  problem,  and  a  detailed  analysis  showing  how  this  can 
be  done  can  be  found  in  [17],  and  [14],  section  2.6.  However,  to  simplify  our 
analysis,  it  will  be  assumed  below  that  the  medium  is  probed  by  a  line  source. 

Then,  the  Fourier  transform  P(x,a))  of  the  pressure  field  satisfies 

[v^+k^n^(x)]  P(x,w)  =  -6(x-i])  (2.1) 

where  k  =  w/c  is  the  wavenumber,  and  n(x)  =  c/c{x)  is  the  refraction  index  of 
the  medium  measured  with  respect  to  some  constant  velocity  c.  Throughout  this 
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paper,  it  will  be  assumed  that  n(x)  does  not  deviate  significantly  from  a 
background  index  nQ{x) ,  so  that 

n^(x)  =  OqCx)  +  U(x)  .  (2.2) 

where  the  scattering  potential  U(x)  is  small.  In  addition,  we  assvime  that 
nQ(*)  is  a  smooth  function  and  that  U(*)  has  a  bounded  support  V,  which  is 
completely  located  on  the  same  side  of  F,  as  shown  in  Fig.  1.  Substituting 
(2.2)  inside  (2.1),  we  can  rewrite  equation  (2.1)  as 

[v^  +  k^nQ(x)]  P(x,u)  =  -  5(x-t2)  -  k^(x)  P(x,«)  (2.3) 

2  2  2 

where  the  operator  Dq  =  v  +  k  nQ(x)  appearing  on  the  left-hand  side  of  (2.3) 
is  the  background  Helmoltz  operator,  and  where  the  second  term  on  the 
right-hand  side  can  be  viewed  as  a  perturbation.  The  solution  of  (2.3)  is 
given  by 


P(x,«)  =  Pq(x,w) 


+  k 


dx'U(x' )P(x‘ ,w)G„(x,x’ ,«) ,  (2.4) 

j  _  _  u  - 


where  the  incident  field  Pq(x,(i))  satisfies  the  unperturbed  equation 

=  -5(x-3)  (2.5) 


and  where  Gq(x,x’ ,w)  is  the  Green’s  function  associated  to  Dq,  i.e. 


DqGq(x,x' ,<j)  =  -6(x-x') 


(2.6) 


By  comparing  (2.5)  and  (2.6)  we  see  that  for  the  experiment  geometry 
considered  here,  the  incident  field  Pq(x,w)  =  Gq(x,i],(i))  . 

The  main  feature  of  equation  (2.4)  is  that  it  is  exact,  i.e.  no 
approximations  are  involved  up  to  this  point.  This  equation  is  known  as  the 
Lippmann-Schwinger  equation  [18],  and  it  puts  in  evidence  the  nonlinear 
relation  existing  between  the  potential  U(x)  and  the  pressure  field  P(x,w). 

In  this  paper,  we  will  linearize  this  equation  by  using  the  Born 
approximation,  whereby  the  total  field  P(x,(o)  is  approximated  by  the  incident 
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field  Pq(x,  u)  =  Gq(x.  t},  w)  inside  the  integral  appearing  in  equation  (2.4). 
In  this  case,  the  scattered  field  P^{x,u)  =  -  Pq{x.w)  can  be  expressed 
as 


Pg(x,{o)  =  J  dx'U{x’)GQ{x.x’.w)GQ(x'.i3.w)  .  (2.7) 

Note  that  the  Born  approximation  is  valid  only  if  the  scattering  potential 
U(*)  is  small  both  in  magnitude  and  spatial  extent  {[19],  Chapter  9).  Then, 
the  relation  between  P^(»,u)  and  U(»)  becomes  1 inear .  and  the  inverse  problem 
that  we  shall  consider  in  this  paper  consists  in  reconstructing  U(x)  for  x  e  V 
from  the  observed  scattered  field  P^Cf.w)  where  ^  =  ]^(s)  is  located  along  the 
array  F.  Note  that  the  receiver  array  F  nay  not  provide  a  total  coverage  of 
the  domain  V,  so  that  in  general  U(x)  will  only  be  partially  reconstructed 
from  the  given  observations. 

Throughout  this  paper,  we  shall  use  the  geometrical  optics  approximation 

^iTr/4 

•^0^-’-  “  1/2  ^  (2.8) 


for  points  x  e  V,  and  for  x'  =  tj  or  x'  e  F,  and  where  for  negative  values  of 
1/2  1/2 

k,  k  =  i(-k)  .  This  approximation  is  a  high-frequency  approximation. 

Equivalently,  it  corresponds  to  assuming  that  the  distance  between  the  domain 
V  where  the  scattering  potential  U(»)  is  located  and  the  source  and  receivers 
is  large  when  compared  to  the  wavelengths  used  to  probe  the  medium,  so  that 

k|x-x’|  »  1  (2.9) 

for  X  e  V  and  x’  =  tj  or  x’e  F.  In  (2.8)  f(x,x’)  is  the  phase  or  travel  time 
function  and  satisfies  the  elkonal  equation 


|v^.p(x,x')  I  =  nQ(x), 


(2.10) 
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auid  the  amplitude  a(x,x')  obeys  the  transport  equation 


av  <P  +  2v  a*v  =  0 

X  X  X 


(2.11) 


along  the  ray  linking  points  x  and  x' .  For  the  special  case  when  the 
background  refraction  index  is  constant,  i.e.  when  nQ(x)  =  we  have 
Gq(x,x',{ij)  =  iH^^^{knQ|x-x'  |)/4  where  denotes  the  Hankel  function  of 

order  zero  and  type  one,  and  the  approximation  (2.8)  reduces  to 

1 


iir/4 


GQ(x,x*,a))  —^72 - ; - n72 

k  2(2TmQ|x-x'  I) 


ikn^lx-x' I 


(2.12) 


Then,  substituting  the  approximation  (2.8)  inside  (2.7),  and  setting 
X  =  £,  we  obtain 


where 


w)  = 

dx'U(x'  )A(x’  .f)e*^^- 

(2.13) 

V 

V 

A(x,f)  =  a(x,f)a(x.T3) 

(2.14a) 

$(x,£)  =  <p(x,g)  +  <p(x.r])  . 

(2.14b) 

F(£.k)  =  P^(£.a))/ik 

(2.15) 

Denoting 


and  letting  f(£,r)  be  the  inverse  Fourier  transform  of  F(g,k)  with  respect  to 
k,  we  find  that 


>r/c 


0 


Ps(^«'r)dT 


f(?.r)  =  - 

f  dx'U(x')A(x',£)6(r-$(x’,£)) 
Jv 


(2.16) 


where  the  time  function  Pg(fit)  is  the  scattered  field  observed  at 

The  Identity  (2.16)  indicates  that  f(^,r)  can  be  viewed  as  a  weighted 
projection  of  the  potential  U(*)  along  the  curve  ^(x,^)  =  r,  and  that  f(f,r) 
is  obtained  by  integrating  the  scattered  field  Pg(f.*)-  Thus,  the  inverse 
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scattering  problem  that  we  consider  here  can  be  formulated  as  a  generalized 
Radon  transform  (GRT)  inversion  problem.  This  was  the  point  of  view  adopted 
by  Beylkin,  Miller  and  Oristaglio  [l]-[3],  who  were  then  able  to  obtain  a 
backprojection  operator  to  reconstruct  the  potential  U(»)  partially  from  its 
projections  f(f,r)  where  £  e  F  and  0  <  r  <  “.  However,  one  disadveintage  of 
the  GRT  inversion  method  is  that  to  implement  the  backprojection  operator,  it 
requires  the  computation  of  the  amplitude  a(x,£)  and  phase  >f>(x,^)  for  every 
point  X  in  the  medium  and  every  receiver  f  e  F,  as  well  as  the  computation  of 
a(x,2)  and  <f’(x,T2)  for  all  points  x,  where  the  position  tj  of  the  source  is 
fixed.  This  amount  of  computations  is  very  large,  and  the  inverse  GRT  is 
therefore  very  costly  to  implement. 

The  objective  of  this  paper  is  to  obtain  an  alternative  reconstruction 
method,  whose  computational  requirements  will  be  smaller  for  a  significant 
class  of  problems. 

3.  The  Backpropagated  Field 

The  inversion  method  that  we  propose  relies  on  the  concept  of 
backpropagated  field,  which  was  first  introduced  in  the  context  of  holographic 
imaging  by  Porter  [20].  This  idea  was  subsequently  used  by  Bojarski  [21]  and 
Esmersoy  [14]  to  study  inverse  source  and  inverse  scattering  problems,  and  it 
was  applied  to  the  solution  of  the  constant  background  inversion  problem  for 
an  impulsive  point  source  and  for  a  plane  wave  source  in  [15]  and  [16], 
respectively.  The  migration  methods  [22]  -  [24]  which  are  currently  used  in 
exploration  geophysics  also  rely  on  the  concept  of  extrapolated  field  to  image 
the  main  reflectors  contained  in  a  scattering  medium.  By  migration  we  mean 
here  a  technique  which  is  used  to  image  the  discontinuities  of  the  scattering 
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potential  U(x).  Migration  differs  from  inversion  by  the  fact  that  in 
migration  the  objective  is  only  to  detect  discontinuities  of  U(x),  whereas 
inversion  methods  seek  to  obtain  some  precise  quantitative  information, 
perhaps  only  partial,  about  the  values  of  the  function  U(x)  or  of  its  Fourier 
transform. 

For  the  experiment  geometry  considered  here,  the  backpropagated  field  is 
defined  as 

PgCx.a))  =  ds  W(a,)P^{^,a))  V^GQ^x,£,a,).i;{£)  (3.1) 

where  W{w)  is  a  filter  to  be  determined  later,  and  where  i3{f)  is  the  unit 
vector  perpendicular  to  F  at  point  g,  oriented  in  the  outwards  direction,  as 
shown  in  Fig.  2.  Note  that  in  integral  (3.1),  the  receiver  position  £  =  £(s) 
is  a  function  of  the  arc  length  s,  but  for  simplicity  this  dependence  will  not 
be  indicated  explicitly  in  the  mathematical  expressions  that  we  shall  consider 
below,  except  at  places  where  our  analysis  will  need  to  take  this  dependence 
into  account.  Here,  v^GQ(x,g,(o)*u{g)  denotes  the  Green’s  function  of  a 

doublet,  i.e.  it  can  be  implemented  by  two  impulsive  sources  of  opposite  signs 
located  perpendicularly  across  the  curve  F  at  a  very  small  distance  from  each 
other,  as  shown  in  Fig.  2.  The  fact  that  we  select  the  complex  conjugate  of 
VgGQ(x,j[,(i))  indicates  that  the  waves  propagate  anticausallv  (backwards  in 

time),  since  we  want  to  reconstruct  the  pressure  field  inside  the  medium  at 
earlier  times. 

The  expression  (3.1)  for  the  backpropagated  field  Pg(x,0)  has  the 
following  physical  interpretation:  Pg(x,(o)  is  the  acoustic  field  which  is 
obtained  by  replacing  the  receivers  along  the  curve  F  by  doublet  sources,  and 
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by  exciting  the  doublet  located  at  point  £  &  F  with  the  source  wavelet 

W(cj)  P  (g  ,u).  This  source  wavelet  is  obtained  by  passing  the  scattered  field 

observed  at  g  through  the  filter  W{u).  From  a  practical  point  of  view, 

P^(x,u)  can  be  computed  in  a  number  of  ways.  The  method  that  we  propose 
consists  in  observing  from  (3.1)  that  P^Cx.w)  satisfies 

DqP^(x.(j)  =0  (3.2) 

for  X  €  r.  Then  if  we  specify  an  appropriate  boundary  condition  for  this 
equation,  the  backpropagated  field  Pg(x, t),  where  Pg(x, t)  denotes  the  inverse 
Fourier  transform  of  P^(x,(o),  can  be  computed  backwards  in  time  by  using  a 
finite-difference  scheme  of  the  type  described  in  [25].  Note  that  except  for 
the  introduction  of  the  filter  W(a)),  the  procedure  described  above  for 
computing  p^(x, t)  is  identical  to  the  reverse-time  finite  difference  migration 
method  for  unstacked  data  which  was  proposed  in  [26],  where  the  boxmdary 
condition 

=  PgCS-t)  (3.3) 

for  ^  e  r  and  t  >  0  was  selected. 

The  only  element  that  has  been  left  unspecified  above  is  the  choice  of 
boundary  conditions  for  p^(x,t).  Strictly  speaking,  to  obtain  an  exact 
boundary  condition  for  P  (x,(i))  one  would  need  to  use  the  expression  (3.1)  to 

/X# 

specify  the  value  of  Pg(x,t)  over  a  boundary  F  located  close  to  the  array  F. 
However,  it  is  more  convenient  to  note  that  when  equation  (3.1)  for  P^(x,cj)  is 
equated  to  expression  (3.11)  below,  and  when  F  surrounds  the  domain  V,  it  was 
shown  by  Esmersoy  in  [14],  pp.  99-105  that 

Pe(£-')  =  Ps(£'‘) 


(3.4a) 
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for  £  fc  r  and 


t  >  max  (<p(x,T7)  -  <p(x.£))  .  (3.4b) 

xfcV 

f 

where  p  (£.t)  denotes  the  inverse  Fourier  transform  of  the  filtered  scattered 
s 

field  W(w)P^(f.w). 

The  Identity  (3.4a)  specifies  a  boundary  condition  for  the  extrapolated 
field  Pg(x, t)  over  the  time  window  (3.4b).  Since  the  inversion  method  that  we 
propose  in  section  4  requires  the  knowledge  of  P^(x,  t)  only  at  t  =  <p(x,i7), 
where  <f>(x.7j)  is  the  travel  time  from  the  source  tj  to  point  x  e  V,  the  time 
window  (3.4b)  is  in  general  sufficient  to  compute  the  values  of  P^(x, t)  that 
we  need.  In  addition,  the  boundary  condition  (3.4a)  is  so  simple  that  even 
when  r  does  not  surround  the  domain  V,  it  may  be  worth  using  it  to  compute  an 
approximation  to  Pg{x. t).  This  scheme  has  given  good  results  in  [14],  [16]. 
However,  as  indicated  above,  a  more  rigorous  method  would  need  to  start  from 
the  definition  (3.1)  of  P^(x,w). 

A  representation  of  the  backpropagated  field  which  will  be  useful  in 
subsequent  derivations  can  be  obtained  by  noting  that  within  the  geometrical 
optics  approximation 

iTr/4  . 

v^Go(x,g,Q)  =  f -  (Vga(x,f)  +  ika(x,£)  v^<p(x,£))e^^’^^^’^^  (3.5) 

kl/2 

Then,  if  we  make  the  additional  approximation 

|Vga(x,f)|  «  |ka(x,g)|  (3.6) 

which  is  of  the  same  nature  as  the  geometrical  optics  approximation,  we  find 


that 


(3.7) 
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The  vector  v(x,g)  =  Vg<p(x,g)  is  tangent  to  the  ray  liiiklng  x  and  g.  as  shown 
in  Fig.  2,  and 

|v(x,£) 1  =  n^Cf),  (3.8) 

so  that 

-  -  k^^^e“^^^\Q(f)cos/3(x.£)a(x.£)e^^'^^-’^^ 

=  iknQ{£)cosP(x.f)GQ(x.£.<u)  .  (3.9) 

^  y\ 

where  p(x,£)  is  the  angle  between  the  normal  vector  i)(£)  and  v(x,£). 
Substituting  (3.9)  inside  (3.1),  and  taking  into  account  the  definition  (2.15) 
of  F(x,k).  this  gives  the  following  expression  for  the  extrapolated  field 

PJx.a,)  =  k^(w)  ds  F(£.k)nQ(£)cosj3(x.£)GQ’‘(x,£.cj)  (3.10) 

which  will  be  used  in  next  section. 

It  is  worth  noting  that  the  definition  (3.1)  of  the  backpropagated  field 
differs  from  the  definition  selected  in  [14]  and  [16],  where  P^(x,w)  was  given 
by  the  Kirchhoff  integral 

P^(x,a))  =  W((o)  ds[P^(£,o))VgGQ^(x,g,<u) 

-  •  (3.11) 

The  motivation  for  considering  the  expression  (3.1)  instead  of  (3.8)  is  that 
it  is  simpler,  so  that  our  subsequent  derivations  will  be  easier  to  follow. 
There  is  however  no  practical  reason  why  (3.1)  should  be  used  instead  of 
(3.11).  An  important  difference  between  (3.1)  and  (3.11)  is  that  the 
definition  (3.1)  requires  only  the  knowledge  of  the  scattered  field  Pg(f ,(o) 
along  the  receiver  array  F,  whereas  (3.11)  requires  also  the  knowledge  of  the 
normal  derivative  ^  Pg(^.“)  along  F.  Since  this  derivative  cannot  usually  be 


measured,  it  may  appear  at  first  sight  that  the  definition  (3.11)  cannot  be 

implemented  directly.  However,  as  was  already  observed  above,  it  was  shown  by 

Esmersoy  [14]  that  with  the  definition  (3.11),  Pg(x, t)  can  be  computed  in  the 

time  domain  by  using  the  homogeneous  equation  (3.2)  with  the  boundary 

condition  (3.4),  so  that  the  normal  derivative  9  P  (^,w)  is  not  needed. 

an  ® 

An  argument  which  indicates  that  the  definitions  (3.1)  and  (3.11)  are 

approximately  equivalent  is  as  follows.  First,  substitute  the  approximations 

r  ik<f  (g) 

P^(£,m)  =  k  \(£)e  ®  (3.12a) 

r 

GQ(x,f,co)  =  k  ®a(x,g)e^^'^^-’^^  (3.12b) 

inside  (3.11),  and  take  into  account  (3.9),  so  that 

r  +1 

a  =  ik  ®  cosP^(g)P^(g,a))  (3.13a) 

an 

5  GQ(x,f,{o)  =  ik  ^0^^^  ‘^osa(x,g)GQ(x,g,w)  (3.13b) 

an 

where  /3  {g)  and  P(x,g)  denote  respectively  the  angles  that  Vj.<p  (f)  and 

Vg>f)(x,g)  make  with  the  normal  to  F.  The  integrals  along  F  of  the  first  and 

second  term  in  (3.11)  are  now  expressed  as 

Ej  =  ds  nQ(£)a^(g)a^(x,g)cosa(x,£)e^^^'^s^^^  ‘f(5>£))  (3.i4a) 

^2  ^  "ik^J^ds  nQ(g)a^(g)a^(x.g)cosP^(g)e^^^‘^s^^^~'^^-'^^^  (3.14b) 

where  r  =  r^  +  r^  +  1.  The  only  difference  between  these  two  integrals  is  the 
weighting  factors  cos/3(x,g)  and  cosP^{g)  which  appear  in  (3.14a)  and  (3.14b) 
respectively.  If  the  method  of  stationary  phase  is  used  to  evaluate  these 
integrals,  we  find  that  the  stationary  points  are  given  by 
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where  t(g)  denotes  the  unit  vector  tangent  to  the  array  F  at  point  But 


and  v^<p{x,£)  satisfy  the  eikonal  equation  and  have  therefore  the  same 

magnitude  nQ(^).  In  addition,  they  are  both  oriented  in  the  outward  direction 
with  respect  to  F,  so  that  at  a  stationary  point  we  must  have 

=  Vg<p(x.£)  (3.16) 

which  in  turn  implies  that 


cos/3g{f)  =  cos/3{x,£)  (3.17) 

Thus,  at  stationary  points  of  (3.14a)  and  (3.14b)  the  weighting  functions 
appearing  in  the  two  integrals  are  equal,  so  that  and  Eg  make  equal 
contributions  to  the  leading  order  asymptotic  expauision  of  P  (x,(i)). 
Consequently,  except  for  a  factor  2  which  can  be  incorporated  in  W(a)) ,  the 
expressions  (3.1)  and  (3.11)  are  approximately  the  same.  Note,  however,  that 
our  analysis  is  only  approximate  since  it  relies  on  high  frequency 
asymptotics.  It  turns  out  that  there  exist  several  special  cases  for  which 
the  two  terms  in  (3.11)  are  exactly  equal.  This  is  the  case  for  example  when 
the  background  medium  is  constant  and  the  receiver  F  is  a  straight  line  [19]. 
4.  Inversion  Method 
4. 1  Reconstruction  Technique 

The  inversion  procedure  that  we  propose  is  just  an  extension  of  a  method 
which  was  developed  earlier  by  Esmersoy  [14]  for  solving  the  linearized 
inversion  problem  for  a  point  source  with  a  constant  background.  The  first 
step  in  our  method  is,  for  every  point  x  in  the  medium,  to  image  the 
backpropagated  field  Pg(x,t)  at  the  source  travel  time  t  =  <p(x,2])/c.  Here 
‘P(x,]])/c  represents  the  time  needed  for  waves  originating  from  the  source  q  to 


reach  x.  This  gives  the  image 


0(x)  =  P„{x.  ^  ‘P(x.i2))  =  ^ 


2ir 


f*“  .  <0  ,  X 

d..  P  (x,o,)e-'  ? 


(4.1) 


e'-  c 

The  image  0(x)  can  be  viewed  as  the  migrated  imeige  obtained  by  applying  the 

source  travel-time  imaging  rule  to  the  extrapolated  field  (see  [26]  for  the 

description  of  a  migration  technique  based  on  this  principle). 

Then,  changing  the  integration  variable  from  (o  to  k=cj/c,  taking  into 

account  the  representation  (3.10)  for  P^(x,a))  and  identities  (2.13),  (2.15), 

and  using  the  geometrical  optics  approximation  for  Gq(*, •,&)),  we  obtain 

0(x)  =  I  dx’  N(x,x’)U(x’)  .  (4.2) 

Jv  “  “ 


where 


N(x,x’)  = 


2Tr 


f 


dk 


\ 

W(w)  ds  n„(f)cosP(x.f). 
Jt  “ 


Let  now 


A(x*.£)  a(x.£)e^kt^(5’-£)  . 


^(5)  =  v^0(x).v^.p(x.T2) 


(4.3) 


(4.4) 


be  the  reconstructed  value  of  the  potential  U(x).  Since  the  vector  v^<p(x,i3) 


is  tangent  to  the  ray  linking  the  source  jj  to  point  x,  and  has  magnitude 

value  is  obtained  by  taking  the  gradient  of  the  migrated  image 
0(x)  in  the  direction  of  the  ray  linking  to  x,  and  by  scaling  the  resulting 
expression  with  nQ(x)/a(x,T2) •  The  scaling  by  nQ(x)/a(x,T])  is  needed  to  take 
dispersion  effects  into  account.  Then,  from  (4.2)  we  find  that 

U(x)  =  f  dx’M(x.x*)U(x’)  (4.5) 

Jy 


where 
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Our  objective  in  the  remainder  of  this  section  will  be  to  obtain  a  simple 
representation  for  M(x,x’),  which  will  be  used  to  show  that  the  reconstructed 
potential  U(x)  is  a  good  approximation  of  U(x). 

4.2  Representation  of  Mfx.x’l 

The  starting  point  in  our  derivation  of  a  representation  for  M(x,x’)  is 
the  expression  (4.6),  where  N{x,x’)  is  given  by  (4.3).  To  evaluate  v^N{x,x’), 

we  assume  that 

lv^(cos  /3(x,g)a{x.g))  I  «  |kcos  /3(x,g)a{x,f)v^$(x.£)  1  .  (4.7) 

which,  again,  is  an  approximation  of  the  same  type  as  the  geometrical  optics 
approximation.  This  gives 

v/(x,x-).v^<p(x,T})  ds  cos  P(x,f) 

•  A(x’ ,g)a(x.f)v^$(x,f)»v^<p(x,3)e^^^^^-  (4.8) 

Next,  we  note  that 

=  u(x,g)  +  Uj(x)  ,  (4.9) 

where  the  vectors 

A 

Ui(x)  =  V^<p(x.rj)  (4.10a) 

A 

u{x,f)  =  V^<p(x,g)  (4.10b) 

are  tangent  to  the  rays  linking  t]  to  x  and  g  to  x,  respectively,  as  shown  in 
Fig.  2.  These  vectors  are  such  that 

luj(x)|  =  |u(x,£)|  =  nQ(x)  ,  (4.11) 

yN. 

and  in  the  following  the  angular  arguments  of  Uj(x)  and  u(x,g)  will  be  denoted 
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by  0j(x)  and  0(x,g),  respectively.  Thus 


=  nQ{x)(l  +  cosa(x.f)) 


(4.12) 


where  a(x,£)  =  0(x,£)  -  0.(x)  is  the  angle  between  the  vectors  u(x,g)  and 

/V 

Uj(x).  Then,  by  combining  (4.6),  (4.8)  and  (4.12),  we  obtain 


2  A  A 

cn  (x) 

=  -  2^;^(x.rj)  J 


ito 


J,  ,  5/2  iTr/4„,  .. 
dk  k  e  W((i)) 


J  ds  nQ(£)  cos  i3(x.f)- 


(1  +  cosa(x,g))  A(x'.g)a(x,f)e^^l^^^^'’£)  "  .  (4  js) 


This  expression  can  be  simplified  further  if,  following  Beylkin  [1],  we 
make  the  following  approximation 

A(x’ .^)  A(x.£)  (4.14a) 

^(x'.£)  ^  ^(x.£)  +  V^$(x.g)*(x'-x)  (4.14b) 

for  X,  x'  e  V  and  ^  e  T.  A  rigorous  analysis  of  this  approximation  in  terms 
of  pseudodifferential  operators  was  proposed  by  Beylkin.  In  this  context,  it 
was  shown  that  (4.14)  has  the  effect  of  retaining  the  singular  component  of 
M(x.x’),  which  is  nonzero  only  when  x’  is  in  the  vicinity  of  x,  and  of 
neglecting  the  smooth  components  of  M(x,x’).  When  the  smooth  components  of 
dropped  from  identity  (4.5),  only  the  rapid  space  variations  of 
U(x),  such  as  the  location  and  size  of  discontinuities,  can  be  related  to 
those  of  U(x).  Equivalently,  in  the  Fourier  domain,  only  the  high  wavenumber 
part  of  the  Fourier  transform  of  U(x)  is  related  to  that  of  U(x) .  Another 
interpretation  of  appromixation  (4.14)  which  is  perhaps  more  appealing  from  a 
physical  point  of  view  consists  in  observing  that  when  V  is  in  the  far  field 
of  the  receiver  array  F.  the  distance  |x  -  x’ |  between  points  x,  x’eV  is  small 
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with  respect  to  |x  -  £|  where  ^eF,  and  in  this  case  the  Taylor  series 
expansion  (4.14)  is  justified.  Note  however  that  for  geophysical  surveys,  it 
cannot  always  be  assumed  that  V  is  in  the  field  of  F,  and  in  this  case  the 
justification  of  (4.14)  proposed  by  Beylkin  is  more  appropriate. 

The  end  effect  of  (4.14)  is  that  (4.13)  takes  the  form 
cn^Cx)  c/2  iTr/4  P 

M(x,x')  = - gjj: — J  dk  k  e  W(w)  ds  nQ(£)cos  P(x.f)* 

„  ikv  $(x.g)*(x'-x) 

(1  +  cos  a(x,g))  a  (x,g)e  -  .  (4.15) 

Now,  consider  an  infinitesimal  ray  tube  originating  from  x  and  centered  around 
the  ray  linking  x  to  g  e  F,  as  shown  in  Fig.  3.  If  ^’is  an  arbitrary  point 
along  the  ray  linking  x  to  £,  and  if  do  and  do'  are  the  cross-sections  of  the 
tube  at  ^  and  respectively,  we  have  [27] 

a^(x.£')dCT’  .  (4.16) 

But  do  can  be  expressed  in  terms  of  the  element  ds  of  arc  length  along  F  as 

do  =  ds  cos  /3(x,£)  .  (4.17) 

Furthermore,  if  d9  is  the  angular  span  of  the  tube  at  x.  and  if  is  located 
at  a  very  small  distance  p  from  x,  we  have 

(4.18a) 

da*  =  pd9  (4.18b) 

a^(x.r)-8S^  (4.18c) 

where  to  obtain  (4.18c)  we  have  used  the  asymptotic  form  (2.10)  of  a 
homogeneous  Green’s  function,  and  the  fact  that  in  the  vicinity  of  x,  the 
medium  is  locally  homogeneous.  This  gives 

cos  P(x,£)  ds  =  d9  , 


(4.19) 
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so  that  M(x,x')  can  be  rewritten  as 


M(x.x')  =  - 


cHqCx) 

167r^ 


dk 


k5/2  J  jjg  +  (,Qg  a(x,g)). 


d0  ,  ikv^^(x.£).(x*-x) 

•  di-  - 

Now,  consider  the  change  of  variables 

T{x)  :  (s,k)  - »  2  =  kv  ^(x«£(s)) 


(4.20) 


(4.21) 


and  let  C(x)  be  the  image  of  I  x  IR  under  this  transformation.  Since  g 
depends  linearly  on  k,  the  domain  C(x)  is  a  cone  whose  span  varies  with  x.  It 
was  shown  by  Beylkin  [1]  that  the  Jacobian  of  this  transformation  is  given  by 


J(x,s.k)  =  |k|nQ(x)(l  +  cos  a(x,f(s)))  {x.£{s)) 

Consequently,  if  we  select 


W(«)  =  - 


4e 


-lTr/4 


1/2  ’ 


c  Ik  |k 

the  weighting  function  M(x,x')  can  be  expressed  as 

M(x,x')  =  — ^  [  dg  e^2*(5  5)  _ 


(27r)‘ 


C(x) 


(4.22) 


(4.23) 


(4.24) 


which  is  the  desired  representation  for  M(x,x'). 

From  (4.24),  we  see  that  when  M(x,x')  is  viewed  as  a  function  of  d  =  x'-x 
parametrized  by  x,  it  is  the  Inverse  Fourier  transform  of  the  characteristic 
function 

1  for  g  e  C(x) 

0  otherwise 

This  shows  that  M(x,x’)  is  a  partial  reconstruction  of  an  impulse,  where  the 


X  (x,g)  = 


(4.25) 


cone  C(x)  specifies  the  range  of  available  projections.  In  the  special  case 


when  C(x)  =  IR  ,  which  corresponds  to  an  experiment  geometry  where  the 
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receiver  array  F  surrounds  V.  we  have  M(x.x')  =  6(x-x').  so  that  in  this  case 
we  can  conclude  that 

U(x)=U(x).  (4.26) 

The  identity  (4.26)  seems  to  suggest  that  our  inversion  procedure 
recovers  exactly  the  scattering  potential  U(x).  It  is  worth  remembering  at 
this  point  that  identity  (4.26),  as  well  as  the  representation  (4.24)  for  the 
kernel  M(x,x’)  are  only  valid  within  the  limits  imposed  by  the  approximations 
we  have  made.  These  approximations  are  of  course  the  Born  approximation,  but 
also  the  geometrical  optics  expansion  (2.8)  and  approximation  (4.14).  It 
turns  out  that  these  last  two  approximations  impose  some  very  severe 
restrictions  on  our  interpretation  of  (4.24)  and  (4.26).  As  we  observed 
above,  approximation  (4.14)  is  automatically  satisfied  if  V  is  in  the  far 
field  of  r.  Otherwise,  as  was  shown  by  Beylkin  [1],  the  smooth  components  of 
M(x,x')  are  neglected,  so  that  expression  (4.25)  for  the  Fourier  transform  of 
M(x,x')  is  only  valid  when  the  wavevector  g  is  large,  and  therefore  identity 
(4.26)  should  be  interpreted  in  the  Fourier  don».in  as 

U(e)  =  U(2)  (4.27) 

for  large  g. 

It  will  be  shown  now  that  the  geometrical  optics  expansion  (2.8)  imposes 
similar  restrictions  on  the  validity  of  relations  (4.24)  and  (4.26).  To  see 
this,  note  that  since  k  has  been  assumed  large,  there  exists  some  constant  r 
such  that  |k|>  r.  Then  the  transform  variable  g  defined  by  (4.21)  takes 
values  over  a  set  C^(x)  which  is  the  image  of  I  x  {k: |>r}  under  T(x) .  This 
set  is  obtained  by  subtracting  from  the  cone  C(x)  a  region  in  the  vicinity  of 
g  =  0.  The  set  C^(x)  is  characterized  in  detail  in  section  5.  The  main 
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aspect  of  this  characterization  is  that  although  k  is  bounded  away  from  zero, 
the  magnitude  of  p  =  |g|  of  the  wavevector  g  is  not  necessarily  nonzero,  since 
it  depends  on  the  length  of  the  vector  v^$(x,g)  given  by  (4.9).  In  fact,  when 

the  receiver  g  is  located  on  the  other  side  of  V  with  respect  to  the  source  rj 
which  may  occur  for  a  vertical  seismic  profiling  experiment  where  the  source 
is  on  the  surface  of  the  earth  and  the  receivers  along  a  vertical  borehole, 
the  length  of  V^$(x,g)  is  close  to  zero.  Nevertheless,  the  constraint  ik|>r 

has  the  effect  that  the  Fourier  transform  relation  (4.25)  for  M(x,x’)  is 

restricted  to  g  e  C^(x).  In  the  case  of  complete  receiver  coverage,  this  also 

implies  that  equality  (4.27)  between  the  Fourier  transforms  of  U(x)  and  U(x) 

holds  primarily  for  large  values  of  g. 

Thus,  both  approximations  (4.14)  and  the  high  frequency  asymptotics  that 

we  have  employed  have  the  effect  of  restricting  the  validity  of  our  inversion 

method  to  the  high  wavenumber  region. 

To  interpret  the  filter  W(w)  which  is  used  to  process  the  scattered  field 

P^(^,w),  also  observe  that  since  k  =  u/c,  we  can  rewrite 

W((o)  =  4c^''^  .  (4.28) 

(-1(0) 

so  that  W(<o)  corresponds  to  an  integration  followed  by  a  square-root 
integration. 

4 . 3  Summary 

To  conclude,  let  us  review  our  reconstruction  procedure  step  by  step. 

1)  First,  the  scattered  field  Pg(£.(o)  observed  at  the  receivers  is 
filtered  with  W((o) ,  which  requires  performing  an  integration,  followed  by  a 


square-root  integration. 
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2)  The  receivers  are  replaced  by  doublet  sources,  and  the  filtered  time 
traces  are  used  as  source  wavelets.  Then,  for  this  distribution  of  sources, 
the  backpropagated  field  p^{x,t)  is  computed  by  a  finite-difference  scheme. 

3)  The  backpropagated  field  Pg(x,t)  is  imaged  at  the  source  travel  time 
<p{x,r])/c.  This  gives  the  migrated  image  0(x) . 

ys 

4)  The  reconstructed  potential  U(x)  is  given  by  (4.4),  which  is  obtained 
by  taking  the  gradient  of  0(x)  along  the  ray  linking  the  source  -q  to  point  x, 
and  by  scaling  the  resulting  expression  with  nQ(x)/a{x,g) . 

The  inversion  procedure  described  above  requires  the  computation  of  the 
extrapolated  field  p  (x,t)  and  the  evaluation  of  the  phase  <p(x,i7)  and 
amplitude  a(x,i3)  for  all  points  in  the  medium.  Since  the  location  rj  of  the 
source  is  fixed,  this  scheme  requires  the  use  of  a  ray  tracing  algorithm  only 
for  the  source  rj.  By  comparison,  the  inverse  GRT  [l]-[3]  requires  the  use  of 
a  ray  tracing  procedure  not  only  for  rj,  but  also  for  every  receiver  f  located 
along  r.  Our  method  can  therefore  be  viewed  as  having  replaced  the  use  of  a 
ray  tracing  scheme  for  receivers  along  F  by  the  computation  of  Pg{x,t),  which 
can  be  done  in  one  batch  operation.  Instead  of  receiver  by  receiver.  The 
advantages  and  disadvantages  of  our  inversion  method  with  respect  to  the 
inverse  GRT  are  primarily  those  of  finite  difference  schemes  with  respect  to 
ray  tracing  methods.  Thus,  when  the  number  of  receivers  is  large,  and  when 
the  finite-difference  scheme  which  is  used  to  compute  p^(x,t)  does  not  require 
too  many  grid  points,  our  reconstruction  technique  is  likely  to  be  faster  than 
the  inverse  GRT. 

5.  Characterization  of  the  Cone  C(x) 

The  accuracy  of  the  reconstruction  method  which  has  been  obtained  in  the 
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previous  section  depends  in  a  crucial  way  on  the  angular  aperture  of  the  cone 
C{x)  appearing  the  representation  (4.24)  of  the  weighting  fvinction  M(x,x'). 
Indeed,  as  the  angular  aperture  of  C(x)  increases,  within  the  approximations 
we  have  made  M(x,x')  gets  closer  to  an  impulse  6(x-x’),  and  the  Fourier 
transform  U(g)  of  the  reconstructed  potential  becomes  a  better  approximation 
of  the  U(g)  in  the  high  wavenumber  region.  We  will  show  now  that  the  angular 
range  of  C(x)  is  purely  a  factor  of  the  experiment  geometry,  and  does  not 
depend  on  the  reconstruction  technique  that  we  have  employed.  In  fact,  as  was 
already  observed  above,  the  same  cone  C(x)  appears  also  in  the  inverse  GRT. 
Specifically,  it  will  be  shown  that  the  angular  rar^e  of  C(x)  depends  on  the 
degree  of  coverage  of  the  domain  V  which  is  provided  by  the  receiver  array  F, 
and  on  the  relative  location  of  the  source  tj  with  respect  to  point  x. 

The  first  step  is  to  note  from  (4.21)  that  v^$(x,£)  indicates  the 

direction  of  a  ray  inside  the  cone  C(x) ,  so  that  as  g  varies  along  F,  v^^(x,g) 

spans  the  cone  C(x) ,  as  shown  in  Fig.  4.  The  relation  (4.9)  also  shows  that 
is  sum  of  the  two  vectors  u(x,£)  and  Uj(x),  which  have  the  same 

length,  and  have  angular  arguments  0(x,£)  and  9.(x),  respectively.  Since 

ys  y>w 

u(x,g)  and  Uj(x)  have  the  same  length,  the  angle  separating  these  two  vectors 
is  bisected  by  V  0(x,g),  as  indicated  in  Fig.  5,  so  that  the  angular  argument 

of  is  given  by 

'/'(x.f)  =  I  (e(x.£)  +  0^(x))  .  (5.1) 

yv 

Furthermore,  since  the  source  tj  is  fixed,  u^(x)  and  therefore  0j(x)  are  fixed, 
so  that  in  (5.1)  only  0(x,g)  varies  as  f  moves  along  F.  Consequently,  if  we 
assume  that  0(x,^)  varies  over  the  angular  range  »(x)  =  [0j(x),  02(x)]  as  f 
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moves  along  F,  the  angular  span  of  the  cone  C(x)  is 

HlS)  =  [|  (0i(5)  +  0i(5)).  I  +  0i(5))]  •  (5.2) 

I 

An  interesting  feature  of  this  result  is  that  the  total  aperture  of  C(x)  is  ^ 
{02{x)-0j(x)) ,  which  is  half  the  range  of  8(x) . 

To  interpret  the  above  result,  consider  the  special  case  when  the 
background  medium  is  constant.  In  this  case,  it  is  possible  to  determine  more 
precisely  the  effect  of  the  receiver  array  F  on  the  span  ^'(x)  of  C(x) .  To  do 
so,  we  will  use  for  the  array  F  a  model  of  the  type  considered  by  Porter  [20] 
and  Esmersoy  and  Levy  [16],  where  it  is  assumed  that  F  is  asymptotic  to  two 
lines  wiht  angles  and  Ug  with  respect  to  the  horizontal  axis,  as  shown  in 
Fig.  6.  In  exploration  geophysics,  this  model  covers  the  case  when  the 
scattered  field  is  measured  by  receivers  located  on  the  surface  of  the  earth 


(ttj  =0,  02  =  it) ,  or  the  offset  vertical  seismic  profiling  geometry,  where  the 
receivers  are  located  on  the  surface  of  the  earth  and  along  a  vertical 
borehole,  say  to  the  right  of  domain  V  (a^  =  case  when 

the  receivers  are  located  above  V  Hind  along  two  vertical  boreholes  on  both 


sides  of  V  (a^ 


Then,  since  the  background  medium  is  constant,  the  rays  linking  a  point  x 


in  the  medivim  to  the  source  13  and  to  receivers  along  F  are  straight  lines. 
Thus,  the  angle  0j(x)  is  the  angle  of  the  line  from  tj  to  point  x,  and  by 
keeping  track  of  the  lines  linking  points  along  F  to  x,  we  see  that  as  £  moves 
along  r,  u(x,g)  sweeps  a  domain  D  illustrated  in  Fig.  6,  whose  angular  range 
is  0  =  ,a^+Tr2‘  This  angular  range  is  independent  of  x  and  is  completely 

parameterized  by  the  angles  and  describing  the  array  F.  Consequently, 
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the  range 

^(x)  =  [|-  (ttg-Tr+e^Cx))  .  I  {aj-Hr+0^(x))]  (5.3) 

of  cone  C(x)  is  purely  a  factor  of  the  experiment  geometry,  and  depends  on  x 
only  through  the  angle  0j(x)  describing  the  relative  location  of  the  source  13 
and  point  x. 

The  last  point  that  needs  to  be  clarified  is  the  effect  of  the  high 
frequency  constraint  |k|>r  on  the  range  of  values  of  the  wavevector  p  under 
the  transformation  (4.21).  As  indicated  at  the  end  of  section  4,  in  this  case 
g  e  Cj.{x)  where  C^(x)  is  obtained  by  subtracting  from  the  cone  C(x)  a  set 
located  in  the  vicinity  of  g  =  0.  To  characterize  this  set,  note  from  (4.21) 
that 

p  =  Igl  =  lk|  |v^$(x,f)|  ,  (5.4) 

and  taking  into  account  the  fact  that  the  vector  v^$(x,f)  is  the  sum  of  the 

A,  /W 

vectors  u(x,^)  and  Uj(x)  which  have  the  same  length  ,  we  find  that 

kx^(x.£)|  =  2nQ(x)  cos((0(x,f)-0^(x))/2)  (5.5) 

Thus,  in  the  direction  '/'(x,g)  =  (6{x,g)+0j(x))/2.  the  wavevector  g  must  be 
such  that 

p  >  2ttiq(k)  cos((0(x,£)-0^(x))/2)  =  2rnQ(x)  cos(>//-0^(x))  .  (5.6) 

But  p  =  2rnQ(x)  cos('/»-0^(x))  is  the  equation  of  a  circle  centered  at  ru^(x) 
and  of  radius  .  This  shows  that  C^{x)  is  obtained  by  subtracting  from 

the  cone  C(x)  the  points  which  are  inside  a  disk  of  radius  rnQ(x)  centered  at 
rUj(x),  or  which  are  inside  the  symmetric  image  of  this  disk  with  respect  to 
the  origin,  as  indicated  in  Fig.  7. 

An  interesting  aspect  of  the  above  result  is  that  when  =  0.(x)+Tr,  which 

^  2 
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corresponds  to  the  case  where  0(x,g)  =  0^(x)+Tr,  or  equivalently 
u(x,£)  =  -Uj(x),  then  the  origin  0  belongs  to  the  domain  C^{x).  This 
indicates  that  some  coverage  in  the  low  wavenumber  region  is  possible  for 
certain  source-receiver  geometries.  Note  that  when  u(x,g)  =  -Uj(x),  the 
receiver  £  is  located  on  the  other  side  of  V  with  respect  to  the  source  rj,  but 
on  the  same  ray.  This  geometry  is  of  a  tomographic  nature,  and  arises  in 
exploration  geophysics  for  vertical  seismic  profiling  or  borehole  to  borehole 
surveys,  where  in  the  last  case  the  source  and  receivers  are  located  in  two 
different  vertical  boreholes.  By  comparison,  note  that  when  ^  =  0^(x),  which 

yN  yN 

corresponds  to  0(x,g)  =  0j{x)  and  u(x,g)  =  -Uj(x),  the  wavenumber  p  must  be 
such  that  p  >  2rnQ(x),  and  g  is  therefore  bounded  away  from  zero.  This 
special  case  corresponds  to  the  so-called  zero-offset  experiment  geometry, 
where  the  source  and  receiver  coincide.  Such  a  geometry  can  be  used  as  a 
model  for  surface  surveys  in  exploration  geophysics,  where  the  source  and 
receivers  are  located  on  the  surface  of  the  earth  at  relatively  small  distance 
from  each  other.  The  above  observations  indicate  that  in  order  to  gain  some 
information  about  the  slow  space  variations  of  the  scattering  potential  U{x) , 
a  tomographic  experiment  geometry  must  be  employed. 
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6.  Conclusion 

In  this  paper,  we  have  obtained  a  new  solution  to  the  variable  background 
Born  inversion  problem  for  the  case  when  the  medium  is  probed  by  a  single 
point  source.  This  solution  generalizes  a  reconstruction  technique  proposed 
earlier  by  Esmersoy  [14]  for  the  constant  background  case.  In  our 
reconstruction  method,  the  backpropagated  field  Pg(x,t)  is  first  computed  by  a 
finite  difference  sheme,  and  is  imaged  at  the  source  travel  times.  This  gives 
the  migrated  image  0(x) ,  and  by  taking  the  gradient  of  this  image  along  rays 
linking  the  source  to  every  point  in  the  medium,  and  scaling  the  resulting 
image  appropriately,  we  obtain  the  reconstructed  potential  U(x).  When  the 
receiver  array  provides  a  total  coverage  of  the  domain  where  the  scattering 
potential  U(x)  is  concentrated,  the  reconstructed  potential  recovers  exactly 
the  rapid  space  variations  of  U(x) ,  or  equivalently  the  high  wavenumber  part 
of  its  Fourier  transform.  In  general,  within  several  approximations  which 
restrict  the  validity  of  our  reconstruction  procedure  to  the  high  wavenumber 
region,  it  is  shown  that  U(x)  corresponds  to  the  convolution  at  point  x  of 
U(*)  with  a  weighting  function  obtained  by  reconstructing  an  impulse  from  its 
projections  over  a  cone  C(x) .  The  angular  range  of  the  cone  C(x)  depends  only 
on  the  experiment  geometry,  and  not  on  the  reconstruction  method  that  we  have 
employed. 

The  most  significant  computational  requirement  of  the  procedure  described 
above  is  the  computation  of  the  extrapolated  field  p^(x,t)  with  a 
finite-difference  scheme.  By  comparison,  the  generalized  Radon  transform 
inversion  method  [l]-[3]  requires  the  use  of  a  ray  tracing  scheme,  and  the 
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evalxiation  of  the  phase  and  amplitude  along  each  ray,  for  every  receiver  along 
the  receiver  array  F.  We  expect  therefore  that,  for  a  wide  class  of  problems, 
our  method  will  be  faster  than  the  inverse  GRT  but  a  detailed  comparison  of 
the  relative  efficiency  of  these  two  techniques  will  be  necessary.  A  number 
of  factors  are  likely  to  influence  this  comparison.  The  first  factor  is  the 
number  of  receivers  along  F:  clearly,  as  this  number  grows,  the  numerical 
complexity  of  the  inverse  GRT  increases  significantly,  whereas  the  number  of 
operations  required  to  compute  the  backpropagated  field  Pg(x,t)  remains 
approximately  the  same.  Another  factor  is  the  complexity  of  the  background 
refraction  index  profile  .  It  is  known  that  ray  tracing  schemes  perform 

very  well  for  simple  media,  but  are  relatively  inaccurate,  and  difficult  to 
use  for  complex  media.  By  comparison,  the  finite-difference  method  is  a  brute 
force  technique  which  is  not  affected  significantly  by  the  complexity  of  the 
background  profile  nQ(*).  Finally,  another  important  factor  in  our  comparison 
will  be  the  relative  location  of  the  domain  V  that  we  want  to  image  with 
respect  to  the  receiver  array  F.  If  this  domain  is  not  too  far  away  from  F, 
the  finite-difference  technique  is  relatively  easy  to  use,  since  it  does  not 
require  the  computation  of  Pg(x,t)  over  a  very  large  region  of  space. 

However,  if  V  is  far  away  from  F.  the  finite-difference  technique  becomes  very 
slow,  whereas  the  ray  tracing  approach  is  not  significantly  affected  by  the 
increase  in  size  of  the  domain  that  we  consider.  The  above  observations  seem 
to  suggest  that  the  inverse  GRT  and  the  method  that  we  propose  here  are  almost 
complementary,  in  the  sense  that  one  method  will  perform  best  on  problems  for 
which  the  other  is  least  suited.  However,  these  conclusions  are  rather 
tentative,  since  as  mentioned  above,  a  detailed  comparison  of  these  two 
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methods  has  yet  to  be  performed.  Note  that  our  reconstruction  technique  gives 
good  results  in  the  constant  background  case  [14]. 

Several  assumptions  have  been  made  throughout  this  paper.  The  first  of 
these  is  that  the  medium  that  we  consider  is  two-dimensional.  The 
reconstruction  procedure  which  has  been  presented  above  can  be  extended  in  a 
straightforward  way  to  3-D  media.  Another  significant  restriction  which  has 
been  imposed  is  that  the  background  refraction  index  nQ(*)  is  a  smooth 
function.  This  assumption  is  somewhat  unrealistic  in  practice,  since  most 
geological  structures  exhibit  significant  discontinuities.  It  would  therefore 
be  of  interest  to  extend  our  reconstruction  technique  to  the  case  when  the 
background  index  is  discontinuous.  The  main  difficulty  in  this  context 

is  that  at  interfaces  where  nQ(»)  is  discontinuous,  the  incident  waves  are 
partly  reflected,  so  that  multiply  reflected  waves  exist,  and  in  this  case  the 
asymptotic  form  (2.8)  of  the  Green’s  function  is  not  valid.  If  the  reflected 
waves  are  neglected,  and  if  only  the  transmission  losses  at  interfaces  are 
taken  into  account,  it  was  shown  in  [1]  how  to  reconstruct  U(*).  However,  it 
would  be  ultimately  desirable  to  obtain  a  reconstruction  technique  that  takes 
into  account  the  reflected  waves  appearing  in  the  backgrovmd  model.  In 
addition,  we  note  that  our  reconstruction  technique  applies  only  to  before 
stack  data,  i.e.  to  data  collected  from  a  single  experiment.  In  actual 
geophysical  surveys,  several  experiments  are  carried  out  for  different  source 
and  receiver  locations.  It  would  therefore  be  useful  to  find  a  way  of 
combining  the  reconstructed  potentials  Uj(x),  1  <  i  <  k  obtained  for  several 
experiments,  where  k  is  the  number  of  experiments.  Note  that  in  this  case  the 
cones  Cj(x),  1  <  i  <  k  at  a  point  x  will  provide  different  angular  coverages 
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for  different  experiment  geometries.  A  simple  averaging  of  the  potentials 
Uj(x)  does  not  give  a  satisfactory  solution  to  this  problem,  since  this 
average  would  have  for  effect  to  weight  too  heavily  the  regions  where  the 
cones  C^(x)  overlap.  Another  way  of  formulating  this  problem  is  to  consider 
the  zero-offset  geometry,  where  the  data  is  given  for  coincident 
source-receiver  pairs  located  along  a  curved  array.  This  data  is  obtained  by 
applying  the  common  depth  point  (CDP)  stacking  process  [23]  to  the  data 
collected  from  a  large  number  of  experiments.  For  a  2  1/2-D  zero-offset 
geometry,  and  when  the  background  refraction  index  is  constant,  an  inversion 
technique  which  relies  on  the  backpropagated  field  was  proposed  by  Esmersoy 
[14].  We  expect  that  this  inversion  technique  can  be  extended  to  the  variable 
background  case  by  using  ideas  similar  to  those  that  have  been  discussed  in 
this  paper. 
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Figure  Captions 

Figure  1:  Scattering  experiment.  The  medium  is  probed  by  an  impulsive  point 
source  located  at  tj  and  the  scattered  field  is  observed  along  the  receiver 
array  F. 

Figure  2-  Parametrization  of  the  rays  linking  point  x  to  the  source  rj  and  to 
receiver  The  doublet  sources  generating  the  backpropated  field  are 
indicated  by  +  and  -  signs. 

Figure  3:  Infinitesimal  ray  tube  originating  from  point  x  and  centered  around 
the  ray  linking  x  to  receiver  g. 

Figure  4:  Cone  C(x)  spanned  by  v^$(x,£)  as  £  moves  along  F. 

Figure  5:  Construction  of  v^^>(x,g^)  by  summation  of  the  vectors  u(x,£)  and 
u.(x). 

Figure  6:  Receiver  array  with  angular  aperture  [a^,  D  denotes  the 

domain  swept  by  the  vector  u{x,^)  as  g  moves  along  F. 

Figure  1-  Range  C^{x)  of  the  wavevector  g  under  the  high  frequency  constraint 
|k|>r. 
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